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SUMMARY 


The flow over a helicopter rotor blade In forward flight Is an important example 
of three-dimensional time-dependent flow. The boundary layers on the rotor blade set 
loss levels and control retreating blade stall. As a consequence there is considerable 
Interest in developing a numerical scheme for solving the time-dependent viscous com- 
pressible three-dimensional flow equations to aid in the design of helicopter rotors. 

In the present report, the development of a computer code to solve a three- 
dimensional unsteady approximate form of the Navier-Stokes equations employing a 
Linearized Block Implicit technique in conjunction with a QR operator scheme is 
described. Results of calculations of several Cartesian test cases are presented. 

These results indicate that the computer code can be applied to more complex flow 
fields such as these encountered on rotating airfoils. 
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INTRODUCTION 


The behavior of boundary layers on wings and bodies has long been of Interest to 
aerodynamlclsts . In both steady and unsteady flows the boundary layers are known to 
govern a major portion of the losses and to significantly Influence the vehicle lift 
and moment coefficients. When the flow is steady, boundary layer prediction schemes 
based on numerical solution to the governing partial differential equations of motion 
have reached a high level of sophistication and predictive accuracy, even in three 
space dimensions. In unsteady flows, such as are commonly encountered in rotary 
winged aircraft, some progress has been made in two space dimensions but little to date 
has appeared on unsteady three-dimensional boundary layers. 

Two particular problems arise with time-dependent three-dimensional boundary 
layers relative to the steady case. The first of these is the rather obvious one of 
time integration with its added requirements of transient accuracy coupled with an 
increase in the computational labor, The second of these is the so-called negative 
cross flow problem, which to some extent has troubled the steady boundary layer 
prediction schemes. Kendall, et al (Ref. 1) discuss the negative cross flow problem 
for steady three-dimensional be ,ndary layers in a very illuminating fashion. This 
particular problem arises when the spanwise component of velocity changes sign and will 
be discussed in detail subsequently. Because of the interest by external aerodynamicists 
in swept wing boundary layers where the negative cross flow problem (in this case flow 
from tip to root) is not usually encountered, the negative cross flow problem has not 
received a great deal of attention to date. However in transient flows, particularly 
those encountered on rotor blades in forward flight, negative cross flows are frequently 
encountered. For instance, the advancing rotor blade has cross flows of one sign during 
the first ninety degrees of rotation and these can change sign over part of the blade 
during the second ninety degrees. 

Thus to be of practical value, time-dependent three-dimensional boundary layer 
prediction schemes require high computational efficiency and transient accuracy 
coupled to the ability to treat arbitrary cross flow profiles. 

In this report we describe the development of a computer code for the efficient 
solution of three-dimensional time-dependent viscous flows on fixed and rotary aircraft. 
The Linearized Block Implicit (LBI) technique of Briley and McDonald (Ref. 2) in 
coordination with a tridiagonal QR operator scheme (Ref. 3) is employed to solve the 
reduced turbulent Navier-Stokes equations which are derived for nonorthogonal coordinates 
in generalized tensor form. The rationale for the choice of this approach is discussed 
in detail in Ref. 3. 


The basic assumptions made in the derivation of these equations are that the 
pressure does not vary normal to the shear layer and that in the energy equation 
the square of the normal velocity is neglected with respect to the other velocity 
components (To = constant) . The latter assumption is included only for computational 
simplification purposes and is not essential in the analysis. A novel method is 
employed for solving the continuity equation in conjunction with the reduced 
Navier-Stokes equations. The continuity equation is split by employing the Douglas- 
Gunn procedure to obtain a consistent approximation to the full equation which is 
then solved as an integral. Results of computations on model problems in Cartesian 
coordinates are presented. 
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ANALYSIS 

Background 

In this section the requirements of a three-dimensional unsteady viscous flow 
computer code for flow over airfoils are discussed. 

Three-dimensional boundary layers occur on the wings and fuselages of both 
fixed and rotary wing aircraft. In both types of vehicles, the boundary layers are 
important in setting loss levels and determining useful operating ranges. As is well 
known, boundary layers are sensitive to pressure gradients. In time-dependent flow 
the temporal acceleration terms appear in the momentum equation in a form very similar 
to the conventional imposed pressure gradient and so for qualitative evaluation purposes 
can be regarded as "pseudo" or "auxiliary" pressure gradients. Viewed in this manner 
the temporal acceleration terms can be seen to influence quantities of practical 
importance such as skin friction, displacement thickness and the onset of separation. 

At the range of frequencies typically encountered in rotary wing aircraft aerodynamic 
problems, it is clear, for instance, from the extensive review of McCroskey (Ref. 4), 
that very significant transient boundary layer effects can be observed. 

In examining the flow problems of practical interest such as loss levels or the 
onset of separation it is evident that all three space dimensions must be considered. 

In conventional aircraft the sweep effect is of Interest and inherently three-dimensional. 
In rotary wi-'g aircraft in forward flight clearly very substantial transient changes 
occur in what might be termed the local sweep angle. However, generally speaking, the 
boundary layers remain thin unless catastrophic flow separation occurs or the flow at 
the wing or rotor tip is considered. As a consequence it might be supposed that the 
usual three-dimensional thin boundary sheet approximations (Nash and Patel, Ref. 5) 
could be used to produce a valid set of governing equations. Fortunately some improve- 
ments in thin boundary sheet approximations are possible as a result of having to 
eliminate the negative cross flow problem mentioned earlier. 

The negative cross flow problem is best explained in a somewhat intuitive manner, 
and a good physical description of the problem is given by Kendall, et al (Ref. 1). 

Looking at the suction surface of a conventional swept back wing the boundary layer 
cross flow, w, is usually outward in the z positive direction along the span from 
root to tip. Thus conventional boundary layer integration schemes have developed by 
forward marching the streamwise velocity u in the streamwise x direction and simul- 
taneously marching out along the span in the z positive direction. In view of the 
physics of the problem, the spanwise marching scheme does not normally encounter 
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negative w, i.e., spanwise Inflow. This is very fortunate because it is difficult, 
indeed it could be argued impossible, to structure a physically satisfactory uncondi- 
tionally stable noniterative scheme which permits forward marching in the spanwise 
direction with a negative w cross flow. At least intuitively the problem of negative 
cross flow implies information being transferred upstream against the spanwise marching 
direction. Conventional stability analyses confirm the inability to forward march into 
regions of significant negative w, From experience with attempts to march the two- 
dimensional boundary layer equations into a region of separated flow and its obvious 
relationship to the negative cross flow problem, it is not surprising that spanwise 
marching into a negative cross flow region is not accomplished without special treatment. 
Recently conventional boundary layer developers have been turning to performing an 
implicit spanwise construction to remove the restriction of only positive cross flows 
(Kendall, et al, Ref. 1). Lin and Rubin (Ref. 6) in their predictor-corrector boundary 
region solutions for flow over a yawed cone at moderate incidence also show that allow- 
ing diffusion in the spanwise direction not only eliminates the problems associated 
with negative cross flow, but improves upon the solutions obtained by three-dimensional 
boundary layer techniques. 

Boundary conditions applied at the tip can influence the flow inboard, if 
required by the physics of the flow. For these reasons the implicit spanwise con- 
struction has been a feature of the three-dimensional duct flow analysis of Briley 
(Ref. 7) and McDonald and Briley (Ref. 8). As a consequence of these observations 
and the need to remove the negative cross flow restriction, a spanwise implicit formu- 
lation seems mandatory for the rotary wing applications and at least desirable for 
fixed wing applications, especially as it can be had for a very modest increase in 
code computational labor. Based on the experience in Refs. 7 and 8, the spanwise 
implicit sweep would only result in about a 20% increase relative to the explicit 
spanwise marching approach. The extension of the conventional three-dimensional 
boundary layer equations to allow spanwise diffusion is easily accomplished, and in 
view of the Improved physical representation which thus follows, it is recommended 
and has been implemented in this effort, 

As a matter of course it has been assumed that normal to the wall an implicit 
formulation would be structured. In recent years for boundary layer type problems 
there has been little dispute as to the efficiency gains to be had from an implicit 
formulation normal to the wall (Ref. 9). However in the streamwise direction for 
steady 2-D flow, the equations are normally forward marched and the implicit stability 
obtained entirely from being implicit in the normal to the wall direction. In time- 
dependent flows a similar structure is to be had so that at each time level one 
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streamwiBe (explicit) forward inarching sweep could be made with two implicit sweeps 
in the spanwise and normal directions to give the desired unconditional «tabllity. 

As mentioned earlier the explicit sweep would probably require less computational 
effort by about 20% than an implicit streamwlse sweep and of course less storage, 

However, since the solution is being time marched the opportunity to take a stream- 
wise implicit sweep at roughly the some cost as the explicit sweep does arise. If 
one does perform a streamwlse explicit sweep, then the linearization of nonlinear 
terms is performed about the known spatial marching level. If an Implicit streamwlse 
structure is adopted, then full time linearization can be utilized. That is the 
linearization of the nonlinear terms is performed about the known time level. As is 
pointed out in Ref. 8, it is easier to obtain a consistent spatial-temporal order 
accurate linearization by marching in time than in space (in time the nonlinear 
marching derivatives have the form u^ whereas in space marching they have the form 
u^u_j). Further by structuring implicitly in the space marching direction, (small) 
regions of axial reverse flow would be permitted. As a result of these combined 
benefits of linearization and separation, a streamwlse implicit structure is advocated 
and has been implemented in this effort. 

Transient calculations mean that, in essence, a full 3-D spatial Integration is 
carried out at each time step. Thus, spatial accuracy is very Important to minimize 
the spatial grid point density for efficiency since many time steps are contemplated 
in a given cycle. In order to get the most out of a given spatial difference formula, 
the errors from representing nonlinear terms by linear combinations of terms should be 
less than or equal to the spatial discretization errors, If the linearization intro- 
duces a greater error than the spatial differencing, then either a coarser spatial 
mesh could be used, or iteration, or some form of linearization Improvement is 
called for. Iteration across a time step is not recommended since this only reduces 
the linearization error and computationally costs as much as a complete time step. 

Cutting back the time step would be preferable to iterating to preserve the ' lineariza- 
tion error at some acceptable level, since cutting back on the time step would Improve 
both the transient error and the linearization error. This point is clearly demonstrated 
in Ref. 3. To obtain a linearization, which introduces errors of at most the same as 
the spatial difference formulae, a Taylor series expansion about the known time level 
can be performed. This process clearly demands a formal block, i.e., coupled, treat- 
ment of the system of equations. For instance in the streamwise momentum equation 
a typical term is linearized: 

, x n 4-1 n+i n . n n+i n n . A *Z\ 

(uw) « u w +u w -u w +0(At ) 
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n+1 

and clearly one cannot lag w at the old time level n without introducing a first 
order time e -’or in order to get an uncoupled system, l.e., w not appearing in the 
streamwise momentum equation. Thus formal linearization and consideration of the 
resulting errors indicate the coupled system ought to be treated from the accuracy 
point of view, inis is further reinforced when it is realized that block, i.e., 
coupled, systems are not computationally expensive (in a relative sense). 

Additionally a second type of approximation arises unconnected with linearization 
but arising from basic coupling terms in the original equations and if Indeed some 
terms in an equation are time lagged in order to uncouple the equation system and 
these terms are of equal Importance to the terms retained, then again an iterative 
updating is called for in order to achieve stability, accuracy and consistency. 

(This could be termed ad hoc equation uncoupling). Blottner (Ref, 9) has shown that 
many iterations around the ad hoc uncoupled set (>10) are sometimes required in order 
to achieve an overall solution accuracy commensurate with the local difference molecule 
accuracy. The linearization technique is described in Ref. 8, together with its applica- 
tion to block coupled splitting schemes. Schemes of this general type are here termed 
"split linearized block implicit" or split LBI schemes, and are reviewed in detail by 
Briley and McDonald (Ref, 2). 

As a general observation, care is required to obtain acceptable transient accuracy 
for long time integration with conventional finite difference schemes. A Crank- 
Nlcolson centered time implicit scheme for Instance, although second order in time, 
shows quite a dispersion problem (relative to other schemes) on the simple pure 
convection problem. However, the problem of transient accuracy is significantly 
reduced in the typical boundary layer problem since the time dependency is continuously 
Input through Initial and boundary conditions and relatively the concern is with "short" 
time integrations. The computational problem is more of what the phase lag of the 
wall shear is, relative to the prescribed free stream disturbance, than concern over 
the convection velocity of a wave in a shear after a long propagation time. The 
interest is in forced oscillations with a minimum scale of the boundary layer thickness 
over a few cycles of the motion, just enough to obtain repetition cyclically. It is, 
therefore, expected that a significant dispersion problem will not arise with a con- 
ventional implicit scheme. 


The governing equations that are considered here are the Navier-Stokes equations, 
continuity, energy and the equation of state which are written in generalized tensor 
form for a bo-ik/ oriented coordinate system (boundary layer coordinates) . In accordance 
with the boundary layer assumptions, the normal momentum equation is eliminated and the 
pressure is specified throughout the viscous layer in its stead. For the energy 
equation constant stagnation temperature T Q is assumed. This assumption is a good 
approximation for the flow fields considered, and is thus included here only for 
purposes of simplification, In the analysis that follows, the full energy equation 
could equally well have been used. Employing the equation of state which relates the 
pressure p to the velocity components u and w by an algebraic equation, the problem can 
be reduced to one involving only the three velocity components, u, w and v and three 
equations, the streamwlse and spanwise momentum equations and the continuity equation, 
Hence, we consider a block-three system rather than a block-four system which leads 
to a significant reduction in computer time. If the full energy equation were to be 
considered, a block-four system would result due to the inclusion of the temperature 
as an additional unknown. 


Coordinate System 

Since the goal of *iMs effort is to solve for the flow over airfoils an under- 
standing of tin tyjr* of geometries to be considered is essential to guide the choice of 
the coordinate system and the structure of the computer code. 

Consider a typical finite span swept wing airfoil as shown in Fig. 1, The coordi- 
nate system is not only dependent upon the geometry of the airfoil but also upon the ap- 
proximations that are made to the governing Navier-Stokes equations. As in boundary 
layer theory we also assume that in the approximate form of the Navier-Stokes equations 
the pressure Is constant normal to the shear layer. Inherent in the assumptions is that 
the shear layer is thin. As pointed out by Howarth (Ref. 10) the boundary layer assump- 
tions lead to the conditions that one coordinate direction must be normal to the body 
surface while the other coordinate directions must lie on the body surface. Furthermore, 
the coordinate lines normal to the surface are straight. These conditions uncouple the 
metric data on the surface from that in the normal direction. Hence the metric data for 
the surface coordinates are functions of the surface coordinates alone, while the metric 
data for the normal coordinate direction are functions of that coordinate alone. 

The choice of the surface coordinates is rather arbitrary and is based on considera- 
tions such as the ease of construction or the grid distribution on the wing surface. 
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In the numerical solution of the flow over an airfoil there are many advantages to be 
gained by the judicious choice of coordinates. The most obvious advantage is that the 
physical boundaries of a flow region can us represented by coordinate surfaces. This 
removes the need for fractional cells in general; hence, the complications and loss of 
accuracy associated with a boundary Interpolation are removed. Another advantage Is 
that a uniform numerical method can be used. The solution can then be performed with 
a fixed number of cells In any given direction and with a uniform mesh spacing. 

In Fig. 2 we can see the advantages of a nonorthogonal grid which conforms with 
the boundaries of the swept wing and covers the entire airfoil over a Cartesian grid 
whi . does not. In addition the coordinate transformation can be constructed to contain 
distributions for pnysleal space mesh points. In this context, the uniform mesh of 
computational space is simply mapped into a suitably distributed mesh in physical space. 
The resolution of large solution gradients is the major objective in the selection of a 
coordinate mesh distribution, as in the resolution of an attached boundary layer. 

Another wore subtle example is the resolution cu large gradients in computational co- 
ordinates due to regions of high curvature on the bounding surfaces. When the trans- 
formation contains the mesh point distribution there is no need to construct the appara- 
tus for the discrete approximation of derivatives on a nonuniform mesh. This results in 
a savings in both computer logic and storage. 

Therefore, in this work a coordinate system is chosen that conforms with the 
boundaries of the physical domain i.e., the wing surface which in general will be 
nonorthogonal. In addition, in order to suitably distribute grid points in regions of 
large gradients, provisions are made for analytical grid transformations (Ref. 1J ) in 
each coordinate direction. 

In view of the type of geometries to be considered and the assumptions made to 
obtain the. approximate form of the Navier-Stokes equations a specialized nonorthogonal 
coordinate system is advocated where the metric tensor which has four independent 
components is given by 


9|J = 


9)1 9|2 


9)2 9 2 2 


o o g 33 
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The subscripts 1 and 2 refer to the directions on the surface of the body while sub- 
script 3 refers to the direction normal to the body. Furthermore, the metric data in 
the coordinate directions on the airfoil surface do not vary with the normal direction, 
l,e. the metric data in a 1 - 2 surface above the body are evaluated on the 
the body surface (Ref. 10), Since we will be dealing with nonorthogonal coordinates 
it is advantageous to derive the equations in general nonorthogonal coordinates employ- 
ing generalized tensors. In Appendix B a brief description of tensor notation is given. 
Further details can be found In Refs. 12 and 13, 

An important feature of the analysis to follow is that the governing equations 
which are derived, under the prescribed assumptions, are invariant for any coordinate 
system or any grid transformation (although, of course, the physical approximations are 
-lordinate dependent). The grid transformations are absorbed into the geometrical 
coefficients, leaving the equations unaltered in form. This point has a considerable 
effect on the development of the computer code. Since the only geometric information 
that must be input is the definition of the metric data and their derivatives, it can 
be contained in one subroutine without modifying the remainder of the code. 

Governing Equations 

In view of the ultimate goal of this program, to solve an approximate form of the 
unsteady three-dimensional Navler-Stokes equations on airfoil shapes, the governing 
equations are derived in general nonorthogonal coordinates and are given in generalized 
tensor notation. We will show that this notation aids in the ordering of the various 
terms in the equations and in many respects simplifies the construction of the computer 
code. 

In the following derivation the governing equations are nondimenslonalized as 

i 

follows, x with respect to the characteristic length L, the velocity with respect to 

2 2 

U^, density, pressure and temperature with respect to p , and /c respectively 

and time with respect to L/U^. Viscosity is nondlmensionalized with respect to 

Continuity Equation 

Consider the continuity equation written in vector form so that it is independent 
of coordinate system i.e., 

dp _ 

37- + V • pq = 0 (1) 
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where p is the density and q is the velocity vector. Impressing the velocity vector 
in a covariant basis 


q - 



( 2 ) 


■j 

where u is the i-th contravariant velocity component and e. is the covariant basis 

i 1 

vector in the x direction. The velocity vector may be expressed in a number of dif- 
ferent forms, each with certain attributes. Here for the moment the velocity vector is 
expressed in a covariant basis, for simplicity. Later the velocity components will be 
transformed into physical components for numerical solution. The contravariant basis 
exhibits variation in its components for instance in slug flow, If the coordinates are 
such that the metric varies. For boundary layer flows the physical velocity components 
are roughly aligned with the coordinates and exhibit no variation with the metric per se. 
As such, it is felt that the actual computations are better performed on the physical 
components. 

The divergence of a vector (cf Appendix B) is given by 


v * pq - ^u k = , + ^u 1 r 


ik 


(d/>U* 


(3) 


ki k j. 

where pu |, is the covariant derivative, pu , , is the partial derivative in the x 

K k ^ 

direction, J is the Jacobian and is the Christoffel symbol (cf Appendix B) . 

In Equation 3 two forms of the divergence are presented, one involving the 
Christoffel symbol or curvature term directly and the other the Jacobian. The former 
Is perhaps more restrictive since it requires additional smoothness of the geometrical 
quantities. However, herein we use either form solely from the point of convenience. 

For the continuity equation we use the form involving the Jacobian while in the momentum 
equations the form involving the Christoffel symbols is employed for the evaluation of 
the explicit (lagged) diffusion terras. Thus the form of the continuity equation which 
is used can be expressed as 


+ 7 (J .' ;uk) 'k = 0 


(A) 
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M omentum Equations 

The momentum equations In vector form can be written as 


2>t 


pitt *p[ $ + (<•*>? ] = 


( 5 ) 


where a is the stress tensor. 

In generalized tensor notation Equation 5 becomes 




at 


e 1 = cr 


Ik 


k e ‘ 


( 6 ) 


The stress tensor is defined as 


cr 


ik 


f . 2 , ±_ .Ik 

" “ L p+ 3 Re A J J Re 


(7) 


ik 


where ]j is the viscosity, p is the pressure, A is the velocity divergence and s is 

the strain tensor. The Reynolds number, Re, is defined as p^U^'L/y^. The strain 
tensor is defined as 


e lk = u 1 


m 


_mk , k 

g + u 


im 


m 


( 8 ) 


mk 


where and g™ 1 are the components of the metric tensor. Employing the fact that 
6 kj = g k ^ and substituting the definition of the strain into the stress tensor 
we obtain 


°- ik = - (p + I wr A ) gik+ ^ u ' L gmk+ mr uk I gm ' 


m 


m 


(9) 


Substituting Equation 9 into the momentum equation and employing the relationship 

kj 


g *| = o 

■k 


( 10 ) 
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we obtain for the i-th momentum equation in the e^ direction 



+ u k u 





k 


+ WUI..+ WUI 


(ID 


In Ref, 3 it was pointed out that the QR Operator scheme requires that derivatives 
in any direction operate on only one variable. In the momentum equation this require- 
ment prevents the implicit treatment of certain diffusion terms that arise due to the 
curvature of the body. Although these terms are often treated explicitly anyway the 
use of standard finite difference techniques instead of the QR Operators would give one 
the opportunity to treat these terms Implicitly, if so desired. However, the use of 
the QR Operator scheme requires these terms be treated explicitly, This together with 
the quasi-linear form of the governing equations are the major limitations that arises 
in the treatment of the approximate form of the Navier-Stokes equations considered from 
the use of the QR Operator scheme. In the usual boundary layer approximations these 

explicitly treated terms would not appear in the equations since they are of 
- 1/2 

order 0 (Re ) or smaller, and should, therefore, be of little consequence. In 
principle, the quasi-linear and, for instance; the full conservative form of the dif- 
ferential equations, are equivalent. In discrete form, various formulations of the 
governing equations exhibit different properties (Ref. 29). In the present problem, 
no distinct disadvantage appears to arise from the required use of a quasi-linear form 
of the governing equations. 

The requirement that derivatives in any direction operate only on one variable 
would be more restrictive in the treatment of the pressure gradient term in the full 
Navier-Stokes equations. The linearization of this term introduces derivatives of all 
the velocity components in a given direction. According to the limitations of the QR 
Operator scheme described above, some of these terms must be treated explicitly. Since 
an explicit treatment of these terms could reduce the stability bound of the calculation 
scheme, an alternate procedure should be considered. This would involve the ad- 
dition of an auxiliary equation relating the pressure gradient term to the derivatives 
of the velocity components and would increase the block size of the system. An 
assessment to the efficiency of such a procedure has not been carried out and further 
work in this area would appear to be warranted. 

In the discussion that follows, we will first split the governing equations into 
an explicit part and an implicit part in accordance with the QR Operator requirements. 
Thereafter, we will cast the resulting equations into "standard form", so that the 
equations can be appropriately linearized and treated with the LBI technique. 
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Since miffed partial derivatives are commonly treated explicitly in orthogonal 
coordinate systems, we will do likewise in generalized nonorthogonal coordinates and 
extend this concept to include mixed second covariant derivatives. All other second 
covariant derivatives are retained as implicit. Although such a procedure would 
automatically treat more terras explicitly than one does for orthogonal coordinates, 
it simplifies the bookkeeping requirements in the construction of the computer code, 
and is thus adopted here. Furthermore, by not splitting up the covariant derivative 
for the purpose of making in implicit, but rather retaining it as a unit could prevent 
Instabilities that may arise due to time splitting. Tills occurs when two portions cf 
one term should cancel identically but cannot due to their being split between two sweeps. 

Diffusion Terms 

Consider the term 

) l k (12) 

If j = k the term is retained as implicit, and if j i k then it is lagged. We will 
consider the case j £ k first. Upon expanding the explicit part of the diffusion 
term it becomes 



j' 


i 

Note that the first term on the right-hand side of the equation is in conservation form. * 

t 

Although the implicit equations are treated in quasi-linear form, for the purpose of I 

evaluating the explicit terms the most convenient representation is used. The implicit j 

terms, with j = k become (note there is no sum on j) r 


| j) j =/x(u ! ,j),j + [ (/iu n ), j Tjn + /xu n , j r'jn -/J-u l , n F n jj] + (/x s' jnpu " 1 + T jj (14) 


where 


s 'ini = rm in r' mJ -r'^r-vr'Ki 


( 15 ) 
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T W[“Vmj 


r ") + 


r in] 


( 16 ) 


ra=i+l, n« s 5~2i and no sura on m and n 


Since involves velocity components and derivatives in directions other than the 
1-th direction, the terra is also treated explicitly. 

Hence the total diffusion terms for the i-th momentum equation is given in 
quasl-llnear form as 


+ [o| l g il ( l ii,j+ 2 ^.r'j|) -Z C| k g kk / J.r 1 | lk ]u l ,j 

+ [|°iWJ u) +s/[ E,c i k flV k r l i]u 1 } 


(17) 


+ Z T,i C,*g» + 1 {Z g" k [ / xu'|„]| k+ 2 9 ni [M^ k |n 3 Ik} 

j=l K-l n=l n=l 

k/n n^l 


where 


and 


10 0 “ 
0 I 0 
0 0!. 


r 2 i i i 



2 I 


L I I 2 J 


(18) 


(19) 


and repeated indicies do not indicate summation. The last two summations can be 
combined into one explicit term, so that the diffusion term for the i-th momentum 
equation becomes 
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and there is no summation for repeated indices, Note that the diffusion operator 
terms for a given direction have been cast into standard quasi-linear form i.e., 


au M +bu x +cu + d 


Convective Terms 


The convective term for the i-th equation can be written as 

^u 1 u 1 | = “ m ^ j ] 


( 21 ) 


which equals when expanded 


-° ui,)l |i = |{ pu'u'.j+j (22) 


and the last term is nonzero only when j = 3. The full momentum equation is obtained 
by substituting Equations (20) and (22) into Equation (11) and treating the pressure 
gradient and velocity devergence as explicit terms. Since the pressure is specified and 
impressed upon the viscous layer, its specification replaces the normal momentum equa- 
tion. Thus, the streamwise and spanwise momentum equations are the only two retained. 
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Energy Equation 

The energy equation employed here states that the stagnation temperature Is 
constant throughout 



( 23 ) 


2 

The generalized tensor notation q Is given by 

2 i J 

q z = u u J g 

i 1 

where u and u J are the contravariant velocity components. Incorporating the assump- 
tions made concerning the coordinate system we employ, i.e. 



we obtain ■ 

q 2 = (u 1 ' 2 g,, + 2u'u 2 g (2 + (u 2 ) 2 g 22 +<u 3 ) 2 g 33 


3 2 

Neglecting the term involving (u ) with respect to the other terms, and defining 
physical velocity components, i.e. 


Up = u'h, 


w p = u‘ 


we obtain 


To = T + T (U P Z+W P Z)+ h^ U P W P 


( 24 ) 


This is the form of the energy equation used. 

Equation of State 

The equation of state assumes a perfect gas and is given by 


y-1 

P ~-pl 


( 25 ) 
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Linearizations 

The following analyses assume a set of linear partial differential equations. 
However, the convective part of the momentum equation and the continuity equation are 
nonlinear, containing terras that Involve the product of density and velocity components. 
In order to overcome this difficulty we employ the linearization procedure (described 
in Ref, 8 and reviewed in Appendix A) to linearize the aforementioned terms by Taylor 
series expansion about the known time level solution. 

The density is first eliminated by employing the equations of state and energy, 
and thereafter the resulting terms are linearized, These terms are of the following 
form 


- (/y>e n+/3 + </e ")/ +/3 


p "^e n 

T n 

[<u') n + 

9 12 
h, h 2 

(u 2 )°] (u') n +I ® 

p n ^ n e n 

T" 

V) n + 

9|2 

h|h 2 

(u') n ](u z )" +/3 


2(T n ~ 

T 0 > +' 

7 V \ p " +/3 l 

T n 

\7-'/ P n J 


(26) 


where all velocity components are the contravariant ones, and 0 is always a velocity 
12 3 

component, (u , u , u ) while iji can be either a velocity component or a derivative of 
a velocity component. In the case of a teim containing pu'*' we set i^ n =i|) n+ ^ = 1. 

It is important to note that in the preceeding equations the contravariant 
velocity components are used. However, as noted in Ref. 14 it appears advantageous 
to solve for the physical velocity components. Therefore, when the governing 
equations are subsequently cast into a form amenable to the application of the LBI 
scheme, they are transformed so that the physical velocity components appear. 
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The Turbulence Model 


We treat the set of three-dimensional ensemble averaged turbulent reduced Navier- 
Stokes equations. Ensemble averaging permits the appearance of low frequency (rela- 
tive to the turbulence) time dependent "mean" flow. It is, therefore, necessary to 
specify a turbulence model suitable for this problem. 

The approach taken in the present effort assumes an isotropic Lurbulent viscosity, 
relating the Reynolds' stress tensor to mean flow gradients. 


Reynolds stress = cr Rey 



(27) 


Using Favre averaging (Ref. 15) the governing equations then ore Identical to the 
laminar equations with velocity and density being taken as mean variables and vis- 
cosity being taken as the sum of the molecular viscosity, y, and the turbulent 
viscosity, y^ 


Spatial Difference Approximations 


QR Operator Notati on 

In this section implicit tridiagonal finite difference approximations to the 
first and second derivatives and to the spatial differential operator will be con- 
sidered. The very versatile QR Operator notation will be introduced, which allows 
as special cases a variety of schemes such as standard second order finite dif- 
ferences, first order upwind differences, fourth order operator compact implicit 
(OCI), fourth order generalized OCI and exponential type methods. Since all these 
schemes are of the same form, a single subroutine which defines the difference 
weights is all that is required to identify the method, while leaving the basic 
structure of the program unaltered. Subsequently, the results of numerical experi- 
ments employing some of these schemes will be presented. The rationale for the use 
of the QR approach in the present problem is discussed in detail in Ref. 3. 

The QR formulation allows for ADI methods and permits the treatment of systems 
of coupled equations, i.e., LBI methods. Although variable mesh schemes can be 
employed within the QR framework, it is believed preferable to use analytic trans- 
formations to obtain a uniform computational mesh, hence attention is restricted to 
uniform mesh formulations. 


21 


ORIGINAL PAGE 19 
OP POOR QUALITY 


The general concepts and notation for two point boundary value problems will be 
Introduced and then the methodology entended to more general linear and nonlinear 
parabolic partial differential equations In one dimension. The extension to multi- 
dimensional problems will also be indicated. 

Consider the two point boundary value problem 

L(u) « a(x)u xx + b(x)u x + c(x)u « 7{x) (28) 

with u(0) and u(l) prescribed. Derivative boundary conditions , although not described 

here, can easily be incorporated into the framework of th<? Q-R operation notation. 

Let the domain be discretized so that Xj “ (j-l)h, j ■* 1, 2,,.*, J + 1, and 

u(x,), F (x.)j S.ruu (x.) and h * 1/J is the mesh width. The numbering conven- 
J j x ,1 J j 

tion was choBen here to be compatible with FORTRAN coding. 

Without loss in generality for a(x) t 0, Eq. (29) can be divided by a(x) so that 
we may treat Instead the following equation 

L(u) « U xx + b(x)u x + c(x)u » f(x) (29) 

where 

b(x) - b(x)/a(x), c(x) ■ c(x)/a(x) and f(x) ■ f(x)/a(x) 

The spatial differential operator is Identified as 

L(u) .u xx + b(x)u,+ c(x)u 


Substituting the finite difference approximations to the first and second 
derivatives 



U, 


U j + I " U j~l 
2h 


u x (xp + o(h 2 ) 


(31) 


D +D~ 

h 2 


Ui 


h 2 


s ) ■ u xx (x J , + othZ > 


(32) 
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Into Eq. (29) and rearranging, we obtain 

u(u) ~ s, + b j Fj ■ [pr - |jr] U H + [ e J ■ -p-] u J + [ 'F' + Ih - ] ■ f l 


or 


£ j > j l)j h + [h 2 Cj - 2] Uj + £ I + “g^*] u j+l " h z fj 


(33) 


where Rc^ * hbj is the cell Reynolds number# 

Equation (33) can be generalized by introducing operator format, i.e., 


r I U M + r j Cu J + r )%l - ^(q]f H 4-qJfj +qjf j4| ) 


(34) 


where the superscripts (-) minus, (c) center, and (+) plus indicate the difference 
weight that multiplies the variable evaluated at the (j-1), (1) and (j+1) grid points, 
respectively, and where the r^'s and q^'s for grid point j are functions of h, b.^, 
bj , bj + ^, Cj_^ , Cj and Cj +1 * Comparing Eqs. (33) and (34) we can identify the r^ 's and 
q j ' s , viz. , 


r j - j - Rcj /2 


rj » h 2 Cj - 2 


fj ■ I 4- R C j /2 


qj * 0 

c . 

q) ■ I 

qj ■ 0 


(35) 


We now define the tridiagonal difference operators Q and R 


r[u,] 


r Ju H + r 'Uj + q U,«| 
o[fj] •qjfj-i +q'fj + qJ*j« 


(36) 
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Noting that L(u) » f and substituting Eq, (36) into Eq, (15) we obtain 

r[Uj] « h z Q [ U( U )j] « h 2 o[ f j] (37) 

*“1 

Alternatively by employing the inverse operator Q an expression for L(u),. can be 

J 

obtained 

L(u), .-^O-'RU, (38) 


For standard control finite differences Q ® Q _1 « I, the Identity matrix, 

(the spatial operator is given explicitly lr, terms of and U j+1 ) so that 

nothing was gained in obtaining Eq. (38). However, in general, for highter order 
methods Q is tridiagonal and Q ^ is a full matrix. Hence Eq. (38) gives us a means 
of expressing the spatial operator for a wider class of difference approximations. 

The formalism in Eq. (38) Is also applicable for first and second derivatives appear- 
ing alcne (cf. Ref. 28). It should be pointed out, however, that Eq. (38) Is not the 
most general formulation since the compact Implicit formulas cannot be combined to 
yield a single scalar equation relating the spatial operator to the function values 
(Ref. 28). 


In Refs. 3 and 16 a technique due to Berger, et al is described for constructing 
fourth order tridiagonal methods which possess a monotonlclty property as the cell 
Reynolds number Is Increased, R c *> ”, We will not repeat It here. However, the result- 
ing Q and R coefficients are given In Table II, 


Another family of schemes that can be expressed In Q-R operator notation are 
the so-called exponential methods. The Idea, originally due to Allen (Ref. 17) 
(independently derived by II 'in (Ref. 18) and McDonald (Ref. 19)) and employed by 
Dennis (Ref. 20), is to set the difference weights so that the numerical solution is 
equated to the analytic solution for the locally frozen constant coefficient equation. 
The Q and R coefficients of this exponential scheme Is given in Table III. This method 
is second order accurate for Rc^Otl) and becomes first order accurate as Rc + M where 
the scheme reverts to first order upwind differencing. 

Another exponential scheme which is uniformly second order accurate was developed 
by El-Mistikawy and Werle (Refs. 24 and 25). The, "exponential box scheme" which is 
incorporated in their solution of the boundary layer equations with strong blowing, 
is based on a spatial operator of the form given in Eq . (29). Berger, et al (Ref. 23) 
derived the counterpart for an operator of the form given in Eq, (30), but with c = 0, 
The Q and R coefficients are presented in Table IV. Although this scheme reverts to 
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second offer upwind differences as Rc *► ro , it does not possess a maximum principle 
analogous to the ordinary differential equation it is approximating as does the 
exponential scheme of Allen (Ref. 17). In Table V a second order central difference 
scheme is presented which adds artificial viscosity to the spatial operator when 
| Rc | >2 so that |Rc| never exceeds 2. This scheme is employed in the solution of 
several model problems and will be discussed in greater detail in the section on 
numerical results. 


Application to Coupled Nonlinear Parabolic Equations 


Before considering the LBI technique, we discuss some of the limitations placed 
on the Q-R operator scheme in solving a system of nonlinear parabolic equations. 
Given a system of m nonlinear parabolic equations in m unknowns, 


m 

S 


I 


°r p 


, n+i n. 
< u lj - u ij> 
At 


“ Nj ^ (U| , Ug t • • •• ti m > Xj , Xg , Xj, t ) 


j - 1,2,..., j+l 


where is a quasilinear spatial operator, the Q-R formalism carries 

directly over provided that for any equation only one independent variable is operated 

upon by the differential operator, For example, 


I 

a(u,w,v) u t 


■ + b(u,v,w)u x + c(u,v,w) 


is allowed since x derivatives of u only appear, while 

I 

olU|W ^y U t * U s« + Mu.v.wju, + C(u,v,w) + d(u,v,w)w x 

is not allowed since x derivatives of both u and w appear. The approximate form of 
unsteady Navler-Stokes equations used here, when written in quasi-llnear form, falls 
within the class of allowable differential operators. Thus, for the problem being 
addressed in the present study the OCI schemes are applicable. Note that within the 
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splitting approach, non ullowable terms in the OCX scheme such as dw above, may be 
split off and treated by a special implicit sweep. Provided care is taken and for 
instance the Douglas-Cunn formalism is achieved too, no particular problem arises 
other than the cost of an additional implicit sweep is incurred. 

Thus multidimensional problems and/or more general equation forms can usually 
be accomodated by a splitting procedure, which reduces the differential operator to a 
sequence of one-dimensional problems which have the appropriate allowable form. 

However, as with standard finite differences, to avoid the cost of additional implicit 
sweeps special procedures must be applied to cross derivative terms, e.g., extrapolation 
or explicit treatment. 

Linearized Block Implicit Scheme 
Consider a system of nonlinear partial differential equations 

A$ f = ¥ (39) 


where I is a vector of unknowns and $ is a source term vector which is a function 
of x , x , x and t. Extension to source terms which are functions of 4) are dis- 
cussed in Ref. (8). jb is a three-dimensional nonlinear differential operator and 
of the matrix A acting on the momentum equations is equal to pi where p is the 
density and I the unity matrix. 

Equation (39) may be entered about n+3 time level, i.e. t n+ ^ - (n+8)At = 
nAt-H3At = t n +CAt, and written 

* n+/3 T ~ n+/3 — n+/3 -rnf/3 
A n+P [ $ ^ 


0 - 8 - 1 is a parameter allowing one to center the time step, i.e., 8=0 corresponds 
to a forward difference, 8 = 1/2 to Crank-Nicolson and 8 = 1 to a backward difference. 

Equation (40) can be linearized by Taylor series expansion in time about the n 
time level by the procedure described in Ref. 8 to give a second order linearization 

A n [5 n+I - $ n ]/At = _/"[3> n+ £-'5’ n ] -2> n 5 n + f W1) 


where 1 is the linearized differential operator obtained from^) by Taylor Series 
expansion in time. 
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The difference between the nonlinear operator .2) and the linear operator 2 is 
defined as M n « jJ) n - 2 n . The intermediate level n + 0 is defined as 


$ 


= /3<I> nH +(i-/3)*" 


( 42 ) 


Using these relationships and dropping the vector superbar for convenience a two-level 
hybrid implicit explicit scheme is obtained 


A n ($ n+ '-$ n )/At = P ^ n ($ n + l -$ n ) + + M n 

nuL Q 

The vector \p represents all of the terms in the system of equations which are 
treated explicitly. More about this will be said later, but for the moment note that 
may be approximated to the requisite ord°.r of accuracy by some multilevel linear 
explicit relationship, or approximated by ij/ 1 with a consequent order reduction in 
temporal accuracy. 


The operator 2 is now expressed as sum of convenient, easily invertible sub- 
operators 2 * 2 + 2 2 + . . . . 2 In the usual ADI framework these suboperators 

are associated with a specific coordinate direction. Further it is supposed that these 
suboperators are expressed in the QR notation introduced earlier. Writing V and 
M n< l> n as a single source term the algorithm is written 


s/3 [i, n +i 2 n +^ 3 n ][^ n+L ^ n ]+[^"+^2+^ n 3 ]4> n + S 


n +/3 


(44) 


To solve this system efficiently it is split into a sequence of easily invertible 
operations following a generalization of the procedure of Douglas and Gunn (Ref. 28) 
in its natural extension to systems of partial differential equations. The Douglas- 
Gunn splitting of Eq. (44) can be written as the following three-step procedure 

A n [ $*.$"] / A t =&/,"($*- <J> n ) + [^, n +i 2 n +^ 3 n ]<I> n + S n+/3 

(45) 

A " [$>**_ At -- £*"[*• - <t>"] + ***_ <*."] + n 3 ] cj>" + S 

A " [$*”$"] /At = /9^ l n [<t>* -<t> n ] +/9 -<3> n ] + /9i 3 n [$***- 4> n ] 

+ [^ + ; 2 " + + s" + 0 
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which can be written in the alternative form 

[a" [■*>'-<£"] = At +^j ] 4> n + Ats n+ ^ 

[A n -AtM n ][^"-1>"] ■ A n [®*-4> n ] < 46 > 

[A n -At^4'’][4>"*' < -4-"] = 

if the intermediate levels are eliminated, the scheme can be written in the so-called 
factored form 

( A n -/SAti’ l n )A n “ l (A n -/3At4 n )A n ~ l (A n -y9AtA n )($ n+l _ c£» n ) = 

0 (47) 

At(-/, n + j? 2 n +j? 3 n )$ n +AtS n+/3 

At this point it becomes necessary to consider the structure of the operators 
J 1 * J! 2 and 1 ^ . It will be recalled from the one-dimensional scalar problem 
that use of the QR format greatly facilitated the introduction of a wide variety of 
spatial difference formulae. It follows that in the extention to multidimensions 
undertaken here it follows that use of the QR formultion results in the appearance of 
the Inverse operator Q ^ with the sub-blocks of the Jt ■£ ^ > £ 3 operators. In 
order to implement the scheme the inverse operator Q ^ must be cleared. Accordingly, 
the scalar operator Q is generalized to the vector operator (L with (diagonal) sub- 
blocks Q^. In this generalization j = 1,2 apply on the momentum equations and j = 3 
applies on the continuity equation. The i subscript is associated with the coordinate 
directions of the ^ operators. The discretization results in one diagonal sub-block 
for each grid point for each of the three Q^. Each intermediate step of the algorithm 
is now premultiplied by the associated with the 1 1 implicit operator. Writing the 
product operator J! ^ as L^, the inverse operators are thus removed and the scheme 
written, once again dropping the vector superscript for convenience 

[q, A n - At/3L, n ][ <$*- $> n ] = AtL, 4> n + AtQ,[«/ z n +^ 3 n ] <£ n + AtQ,S n+ ^ 

r nr r - (48) 

[Q 2 A n - At /3L 2 n ] [<$**- 3> n ] = Q 2 A n [<£* - $ n ] 

[q 3 A n - ?][****-$?] = Q 3 f[***-*"\ 

^ 2 V = 0 2 “« R 2 $" ^ 3 n$n = V' R 3 

cjjn+i .. ** + q( At 3 ) 
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With the removal of the inverse operator Q - ’ 1 ’, the question of the proper inter- 
mediate level solution boundary conditions can be addressed. As is pointed out by 
Briley and McDonald (Ref. 2 ), the proper intermediate level boundary conditions may be 
derived by running through the intermediate steps in reverse order. Defining a boundary 
condition operator after linearizing the appropriate physical boundary condition by 
Taylor series expansion in time as 

b/W')* g(t,<l>") 

and applying this operator to the algorithm defines the boundary conditions as 

B 3 "0 3 A" [$**-<!>"] = [B 3 n 0 3 A n - At 0B 3 "L 3 n ] [**«*- *"] 

Bi ,"0 2 A"^* -<£"] =[b 2 "0 2 a" -At 0BjV 2 n ][***-$"] <«> 


and note that unless B^L*^ commute (an unlikely event except with Dirchilet boundary 
conditions, where B™ = I) the exact boundary conditions cannot be derived. A number 
of possible strategies are possible at this point aimed at various levels of approxima- 
tion to B^L^. For the present the term AtBB"L^[4> - 4> n ] is neglected. This introduces 
an error of order 0 [At (4) - 4> n ) ] into the solution but note that this error disappears 
at steady state where $*** = 4>** = 4>*.. Neglect of the AtpB^L^ [$-4> n J term is of course 
equivalent to applying the physical boundary conditions on the intermediate level 
variables. 

This completes the general derivation of the algorithm and attention is now 
given to the specific forms of the operators Including the rather special form 
of the component operator for the continuity equation. 

It is worth noting that the operator % or J can be split into any number of compo- 
nents which need not be associated with a particular coordinate direction. As pointed 
out by Douglas and Gunn (Ref. 28), the criterion for identifying sub-operators is that 
the associated matrices be "easily solved" (i.e., narrow-banded). Thus, mixed deriva- 
tives and the complicating terms which might inhibit the use of OCI can be treated 
implicitly within such a framework, although this would increase the number of inter- 
mediate steps and thereby complicate the solution procedure. 

An inspection of Eq. (48) reveals that only the linearized operators l", L.” and 
appear. Indeed, the computer code employs this feature by evaluating these three 
operators before the first sweep, storing them and accessing them as needed in the 
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subsequent three sweeps. In addition, the terms arising from the nonlinear terms are 
immediately absorbed into S as they appear, allowing for an efficient evaluation of 
the terms in the differential equation 

The operators 2 » ^ can k e represented in standard form at each grid point, 

l.e. 


/ n d> n s a n (I) 

M f^l Q |l *1. 


+ a n <t> 

12 1,1 


+ a 


13 


$ + a n <E> 

I 14 2 


+ a $ 

IB 3 


(50) 


In Eq. (50) the subscripts of $ indicate the velocity component (associated with the 
corresponding direction and " , " indicates a derivative. The subscripts of the a”j 
refer to the direction (i) and the terra in the equation (j) respectively. Note that 
the equation xs in quasi-linear form, a necessity of the QR operator technique employed 
here since the coefficients of the derivative operators need to be identified. Alter- 
nate schemes have been proposed by Leventhal (Ref. 33) for equations in conservation 
form but are not considered here. In the following section we will describe how this 
entire operator is discretized by employing the QR operator format, and show how 
the discretization is incorporated into the LSI framework in order to solve the system of 
equations (48). 

We consider first, however, the continuity equation. Since it is a first order 
partial differential equation it does not have the standard form of Eq. (49). Further- 
more, p has been linearized and eliminated in favor of the u^ velocity components so 
that the continuity equation has become an equation for the three velocity components, 
and not density. 

As pointed out by McDonald and Briley (8) skillful partioning of the resulting 
matrix can lead to significant decreases in computation time. An Inspection of the 
system of equations under consideration reveals that substantial savings can be 
realized if the equations are partioned appropriately. Due to the use of a boundary 
layer coordinate system, in the two momentum equations the normal velocity appears 
only in conjunction with terms associated with the normal "3" direction. Hence, in 
the first two sweeps, for the streamwise and spanwise momentum equation one is required 
to solve only for the two corresponding velocity components without the need of con- 
sidering the continuity equation. However, on the third sweep where all 3 velocity 
components appear, one must solve all 3 equations. This strategy reduces the solution 
procedure to the inversion of two 2x2 block matrices and one 3x3 block matrix 
rather than three 3x3 block matrices which accounts for the reduction in computation 
time. 


If we were to consider the full Navier-Stokes equations which include a normal 
momentum equation the aforementioned partioning could not be applied since the 
normal velocity would appear in all three sweeps. 
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The question that arises is how to appropriately split the continuity equation, 
since it need only be solved on the third sweep. Here again the Douglas-Gunn formulation 
gives us the answer. The continuity equation written in conservation form becomes, 



i£- + _L 
at j 

IF [ J /> ul ] = 0 

(51) 

linearizing p we obtain in 

increment form 


A n Au n+l + B n A W n+I + 

At/3 d 
J 

■ [v n A n Au n+l + v n B n Aw n *'+/Av nt| ] 


At dr ) 

r 

in + AtP 
J J 

-^ r [(A n + uV)Au n *' + (u n B n )Aw' ,tl ] 

(52) 


+ U0 
J 

jjf [</>" + w n B n ) Aw n+ 1 + ( w n A n ) A u n+l j 



where all the velocity components are the contravariant components u = u , w = u 
3 

and v « u . J is the Jacobian and 

£ U/* v”] 

b" = ■7s-[g 22 w n + g, 2 u n ] 


Approximating equation (52) by employing the Douglas-Gunn procedure as a third 

sweep equation, we obtain a consistent approximation to the continuity equation, i.e. , 

1 2 
the x derivative term is evaluated at the * level and the x derivative term is 

evaluated at the ** level. The values of the intermediate derivative terms are ob- 
tained after the solution of the first two sweeps of the two momentum equations. 

Note that these terms do not contain the normal velocity. We can thus write the 
equation in symbolic form 


i- 

P 

f 

}" 

p 

t 

i 

t 

i- 


i 
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A n Au ntl + B "Aw nt, + -M [ j{aVAu"*' + vVaw^VAv"*'}] 


= s " ~ Z 3 T" ITT [ J { }] -P 


At_ 

J dx 2 


[•»{}] 




(53) 


3 

Since the only term involving v is in the x derivative term, we can directly 

3 

integrate the equation with respect to x , i.e. 


n A "A U " + ' + B"Aw n V 3 + AI^vVau^+vVaw"*' +/Av"*' ] 

/SAt r •,* /3At 


L 


s n - 


[ ] 


[ ] 


(54) 


t/x 


In the next section we describe how this is done very easily via the Q-R operator 
scheme. The concept of integrating directly the continuity equation is not new. 

Davis (Ref. 34) in his coupled procedure for the solution of two-dimensional steady 
boundary layer equation used a trapezoidal rule to integrate the continuity equation. 
Weinberg (Ref. 35 ) (Ref. 36) also used a fourth order Simpson integration scheme to 
solve the compressible boundary layer equations. Such procedures are stable and offer 
a viable alternative tp approximating the derivatives by finite differences. 

Note that conceptually the continuity equation in integrated form is treated on each 
sweep of the Douglas-Gunn splitting, although in actuality this can be viewed as 
having the same form as each sweep and the integration operator can be incorporated 
into the J and difference operators, and as a result the stability and consistency 
of the original splitting is retained. 

Implementation of the LBI Scheme Employing 
the QR Operator Technique 

Consider the third sweep of Eq. (47) in which we solve both momentum equations 
and the continuity equation. The momentum equations are in the form 

[a P - /3Atf n 3 ] A4>*** s A n A$** (55) 
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where A<t>*** is the column vector of unknowns, u, v, w. Here we have implicitly 
assumed that the equtions have been appropriately normalized and that the contra- 
variant velocity components have been suitably transformed into their physical com- 
ponents. Employing physical components, (cf. Ref. 14) leads to a better behaved solution 
since these components are not unduly influenced by geometrical variations. 

For the streamwise momentum equation we obtain 

A$***= Au, 33 + Ug 3 Au, 3 + CI 33 A u + a 43 Aw + a„Av (55a) 


while for the spanwise momentum equation we obtain 


*T\ * 

A«P = Aw, 33 + b 23 Aw , 3 + b 33 Aw + b 43 Av + b 53 Au 


(55b) 


where we have omitted superscript *** from Au, Av and Aw. Now in equation 55a, we can 
approximate 


Au »33 + ° 23 Au »3 + q « Au by 


the operator equivalent, so that 


33 


or' r i 

Ax 3 2 


Au 


(56) 


f n A ,+,**# Q| 'R|Au 
A a<P = — — + a 43 Aw + a 53 Av 


(57) 


Similar approximations are made for Equation (55b) . After substituting Equation (56) 
into Equation (54) and multiplying thru by Q we obtain for the streamwise momentum 
equation 


[Oi/j n - /3 Xr,]au - /SAt Q,a 43 Aw -y9AtQ|a 53 Av = Q,p n Au** 


(58) 


where A = At/Ax 3 


2 
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Similarly for Che spanwlse momentum equation we obtain 


[o 2 /> n - /3XR 2 ] Aw “ /3At<y) 43 Av -/3At<y> M A u = q 2 /Aw** 


( 59 ) 


We employ the same type of procedure for the continuity equation. Since the continuity 
equation involves only first derivatives, so that they can be represented as 

d Q c*Rc 

It? = Ax 3 (60) 


The operators Qc and Rc are constructed to approximate the weights associated with 
either a second order trapezoidal rule or a fourth order Simpson's rule, l.e. 

Trapezoidal rule 



r c = 0 , r c = " 1 • r c + = 1 


Simpson's rule 



The continuity equation thus becomes 


J A n Au + JB n Aw + Q“ , R c [jA n v n Au+ JbVAw + JpAV ] 


RHS 


(61) 
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where RHS contains all the terms due to the linearization procedure and the terms 
evaluated at the * and ** levels. Multiplying thru by Q c and setting m « At/ Ax^ 
we obtain 


[Q c JA n + /3tdR c JA n V n ]Au + [ Q c JB n + ^3wR c JB n V n ] Aw + [ /3u»R c J/a n ] AV = Q C (RHS) 


(62) 


The resulting matrix derived from Equations 57, 58 and 61 becomes a block 3 tridiagonal 
matrix (Q and R are tridiagonal operators) with each sub block taking on the form 


[ Q |jP n -/3XR, ][ --/3 At 0,043 ][-/3AtQ,a 53 ] 


Au 


Q,(Au**) 

[ --/3AtQ 2 b 53 ] [ ®?.P ~ ] [ _ ^AtQ 2 b 43 ] 


Aw 

= 

Q 2 (Aw**) 

[Q c JA n +/3aj R c JA n v n ] [q c J B n + R c J B n v n j [/3tuR c J/3 n ] 


Av 


Q 0 (RHS) 

C 


This matrix is inverted by standard LU decomposition. 
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Boundary Conditions and Initial Conditions 


The type of boundary conditions employed in the solution of the approximate form 
of the Navler-Stokes equations are described in this section. On the body surface no 
slip is prescribed for all the velocity components, At the outer edge of the viscous 
layer the values of streamwise and spanwise velocity components are also prescribed, 
However, the value of the normal velocity component is not set, but rather computed as 
part of the numerical solution as is the practice in standard boundary layer procedures. 

At the inflow boundaries, (upstream) velocity profiles ore fixed, while extra- 
polation conditions are employed at the outflow boundaries (downstream). Further dis- 
cussion of this matter is given in the section of numerical results. 

The intermediate boundary conditions employed on the first two sweeps are the 
physical ones. Since all the multidimensional problems considered here to date have 
been steady, their has not yet been any need for a high temporal accuracy solution, 
and the imposition of physical intermediate boundary conditions did not impair the 
quality of the solutions obtained. These results are in keeping with the analysis of 
McDonald and Briley (Ref. 2) for second order spatial schemes. (The type of schemes 
used to dace in the present study) . Fourth order methods have not yet been applied 
to any of the test problems in this report. 

The question of proper intermediate boundary conditions for fourth order methods 
until recently has not been resolved. Fairweather and Mitchell (Ref. 37) developed 
nonphysical intermediate boundary conditions for a fourth order solution of Laplace's 
equation, and showed that, in general, the use of noncorrected l.e, physical boundary 
conditions leads to a loss in steady state accuracy for their method. As pointed out 
by Fairweather and Mitchell (Ref. 37) their scheme is inconsistent. It is this in- 
consistency that requires one to use appropriately derived intermediate boundary condi- 
tions in order to recover a steady state solution independent of time. However, if a consi 
tent scheme were to be used, e.g. Douglas-Gumt, then physical boundary conditions can 
be applied without any loss in steady state accuracy. A more complete discussion of the 
implementation of ADI schemes with the appropriate intermediate boundary conditions 
for QR operator schemes (including the fourth-order generalized OCI scheme) will appear 
in a forthcoming report. These conclusions generalize the results obtained by Briley 
and McDonald (Ref. 2) for second order finite difference methods to higher order schemes 
and to those schemes that can be cast into a QR operator framework. 
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The Computer Code 


The type of numerical algorithm employed as well as its formulation has a 
marked impact on the structure of the computer code. We need to consider both the 
number of CPU operations aB well as the memory requirements, Usually the number of 
operations can be reduced at the expense of increasing the amount of storage. However, 
for three-dimensional problems the accessible fast (small core) memory becomes a 
severe limitation even without attempting to optimize the operation count. 

The storage requirements for the solution of the approximate form of the three- 
dimensional Navler-Stokes equations for even modest size grids (e.g. 30 x 30 x 30) 
exceed the available small core memory of a machine like the CDC 7600. One must 
then resort to either mass storage devices such as disks or slow access memory 
large core) . 

In using such devices both access time and transfer rates must be considered. 

When small amounts of data are being transferred frequently (what we are considering) 
then access time becomes a significant factor. Therefore, a combination of strategies 
must be employed in order to optimize both access time and transfer rate. As a result 
of these necessary manipulations, the resulting code was far more complex than the 
one that would have been developed specifically for a machine with '’unlimited 1 ' storage 
resources. 

An investigation of the operation count of the LBI scheme in conjunction with 
the QR Operator technique leads to the conclusion that the most significant fraction of 
time is spent in computing the matrix coefficients, i.e, the linearization coefficients 
and difference weights. This amount far exceeds the time required for the matrix 
Inversion. Hence it is worthwhile to optimize the calculation of these coefficients 
and if possible store their values. This procedure was accomplished by storing the 
operator coefficients on J anc * ^ 3 as they were computed in the first sweep on the 
right-hand side of the differential equation. On the second and third sweeps 
and J ” were accessed respectively and were not recomputed. It was for this reason 
that the formulation of the LBI scheme referred to linearized operators J ^'s 
instead ofj&^'s on the right-hand side of the equation. 

The general structure of the computer code will not be described. After the 
input section and the initialization of data e.g. geometry, grid transformations, 
flowfield, etc. the actual construction of the difference operators is begun. The 
first derivatives of the velocity components and viscosity are obtained for the entire 
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flow field and stored for ready access when needed for the computation of the appro- 
priate terms In the governing equations. Thereafter the terms that are to he logged 
(treated explicitly) are evaluated and absorbed into the function S n . 

The operators 2 ^, 2 ^ and 2 ^ are then computed. These are used to evaluate the 
appropriate Q and R coefficients which are then stored for easy retrieval during each 
of the ADI sweeps. 

In the first sweep the matrix resulting from the application of the 2-y operators 
for the streamwise and spanwise momentum equation is solved as a 2 x 2 coupled system. 
The solution of this system, the * level quantities, are then used to construct the 
right-hand side of the second sweep equations and to evaluate the appropriate * level 
term in the continuity equation. At this point the 2 ^ operator is accessed and again 
a 2 x 2 system of equations for the streamwise and spanwise momentum equation is 
solved. The ** level quantities are then used to construct the right-hand side of the 
third sweep equations as well as the appropriate terras in the continuity equation. The 
third sweep equations consist of the two momentum equations and the continuity equa- 
tion, with the 2 ^ operator being accessed from memory. The resulting 3x3 system of 
equations are solved for the three velocity components. 

After the ADI procedure is completed, the thermodynamic quantities, density, 
temperature and viscosity are computed. The procedure is then repeated at the following 
time step. 

It is noteworthy that the scheme just described operates on vectors, I.e. lines 
of data. Therefore, it could show promise for vectorized machines , 


NUMERICAL RESULTS 


In thiB section we describe the numerical results obtained by using the computer 
code described in the previous sections on several Cartesian test problems. The goal 
of these calculations is to validate the computer code. No specific timing experi- 
ments have yet been conducted to evaluate the code's efficiency. As mentioned in the 
previous section, the code was structured for Che calculation of three-dimensional 
problems so that the disk writing and large core data allocation are necessary to run 
the code efficiently in that mode. For two-dimensional problems these artifices, 
although not necessary except for very large mesh problems, are still employed. Hence 
as an observation, it in noted that the code le not as efficient in the two-dimensional 
mode as it could be. 

Since the primary goal in this portion of the effort was to obtain a working 
computer code, the following test cases are considered: 

1. One-dimensional unsteady oscillating flat plate; 

2. Three-dimensional boundary layer without pressure gradient; and 

3. Two-dimensional steady boundary layer with and without pressure gradient. 

These cases check out three basic features of the code, viz., time dependent behavior, 
pressure gradient effects and three-dimensional effects. 

All the calculations employ second order finite difference techniques. Although 
a fourth order generalized OCX option is incorporated into the code (and the code is 
structured particularly to encompass that scheme), to date none of the cases considered 
were run employing this scheme,. In subsequent work fourth order calculations will be 
compared to second order finite difference as they relate to efficiency gains and 
convergence rates. 


Oscillating Flat Plate 

The first case considered is that of an oscillating flat plate in a unif o , . ip stream. 
It is a one-dimensional unsteady problem and is thus only a function of the normal 
distance to the wall and time. 

Schlich^l.ng (Ref. 38) gives the exact solution to this problem, which he terms 
'Stokes second problem’, as 

u(77,t) = Uqop "Vcosiwi -7)) 
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where n * y \f 'ui and w ■ the frequency of oscillation. 

2v 

The boundary conditions are 

at the wall 

Tj - 0 u(0,t) = UcoCOStut 

v(0,t) = 0 

at the outer edge 

77 -*■ CD u (0,t ) ~ 0 

Since we must specify the velocity / a finite value of n we set u at q - r) to the 
exact value of u(n e ,t). As initial conditions, the velocity field is set to the exact 
value at t = 0. Since the solution does not take into account any initial conditions 
of the fluid, the relationship for u(n,t) represents a "steady state solution". 

The computer is run in the two-dimensional mode, with 5 grid points in the 
streamwise direction, so that additional boundary conditions are required at the inflow 
and exit planes. This is accomplished by setting the first derivative of the velocity 
components to zero there, which eliminates any variation in the streamwise direction. 

For the problem considered, the following parameters are employed. 

a ) = 100 rad/sec 

Re = 100 

n - 5.0 
e 

A uniformly distributed grid of 11 points in the normal direction is used as well 
as a fixed time step of o>At « 5°. Second order central differences are employed for 
the calculation which includes all the viscous terms appearing in the governing equation. 
Since there is no variation in the streamwise direction and the normal velocity is zero 
throughout, the governing equations reduce computationally to their one-dimensional 
counterpart. In Fig. 3 the computed skin friction coefficient is compared to the 
analytical result over one and one-half time cycles. The analytical value is 
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Even with a relatively coarse meull there is good agreement between theory and the 
numerical results. Although the Crank-Nicolson results (3 =■ 1/2) are more accurate 
than the fully implicit calculation (3*1) as can be seen in Fig, 4 where the L ^ 
error of each calculation is compared, the skin friction results appear much closer. 
We can account for the discrepancy by looking at the evaluation of the skin friction 
coefficient. 

The error in the skin friction coefficient is made up in part from the error xn 
the solution of the differential equation and in part from the approximation of the 
first derivative in the computation of the shear. The error in the evaluation of the 
first derivative is given as 

err0r (u^so) = <*,U VVV = a 2 (v^W) 


where and 02 are bounded constants 

For our case = 50 so that the total truncation error is dominated by this term. 

The results of this calculation indicate that accurate time dependent calculations 
can be obtained. 


Three-Dimensional Flat Plate Boundary Layer 

The second problem considered is the three-dimensional viscous flow over a flat 
plate skewed at an angle of 45° to the flow direction. In Fig. 5 the computational 
domain is shown in which the solution is obtained. An equally spaced mesh of 11 
points in each of the three directions is employed. 

Two cases are considered. In the first case boundary layer assumptions are 
made, and the only viscous terms retained in the momentum equation are those that 
appear in the standard three-dimensional boundary layer equations. Diffusion terms in 
the streamwise and spanwise directions are ommitted. In addition, first order approxi- 
mations to the streamwise and spanwise first derivatives (in the marching directions) 
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are made. The second case considered retains all the viscous terms. The derivative 
approximation to the spatial operators employed entered differences with artificial 
viscosity (cf Table V). 

The boundary conditions applied are no slip at the surface and velocity specified 
at the outer edge of the viscous layer. There the velocity components are set to 1 fj 2 
Note that the vector sum of the two components is 1. At the inflow boundaries 1 and 2 
(cf Fig. 5) the velocity profiles are set while at the outflow boundaries 3 and 4 
second derivative extrapolation conditions are employed. 

Both cases compared well with one another and to the Blasius solution. In Fig, 6 
the displacement thickness 6* and the momentum thickness f> along the diagonal of the 
computational domain are compared to the theoretical result. Although symmetry was 
retained across the diagonal, at a constant £ value, there was some variation of the 
integrated properties as a function of £ from their values on the diagonal. 

Two-Dimensional Howarth Flow 

The third case considered is a two-dimensional steady flow problem; the laminar 
boundary layer on flat plate in the presence of an adverse pressure gradient, the 
Howarth flow. Here again we consider two cases. The first employs the boundary layer 
approximations while the second retains all the terms in the streamwise momentum 
equations . 

A 21 x 21 mesh is used in both cases. The boundary conditions employed are 
no slip at the wall, streamwise velocity specified at the outer edge, velocity profile 
specified at the inflow plane and second derivative extrapolation at the downstream 
plane. 

The external velocity field is linearly retarded, i.e. u = 1 - Kx, so that 
there is an adverse pressure gradient. Due to the adverse pressure gradient separa- 
tion will occur at some point downstream. Howarth (Ref. 39) computed the separation 
point at a value of x* = Kx = .1200. In the calculations considered here we choose 
our domain to span .05 - x* - .12 which terminates at a location corresponding to 

Howarth' s predicted separation point. 
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The initial conditions for the problem were constructed from the local Falkner 
Skan similarity solution with the appropriate pressure gradient parameter imposed 
at the corresponding streamwise location. 

In Fig. 7 the computed skin friction coefficients for both cases are compared to 
the results of Howarth. The case considered corresponds to Briley’s cases (Ref. 40), 
l.e. Re ■ 62,500 and k « 3. There is good agreement over most of the domain, with an 
overprediction of the shear as the separation point is approached. This discrepency 
is not surprising due to the coarse mesh employed. It is interesting that the 
"boundary layer" case gives somewhat better predictions near the separation point 
than the "full equation" case. This is probably due to the lack of upstream influence 
in the boundary layer solution which did not allow the flow to adjust upstream of 
separation. Both calculations converged within 30 - 40 iterations, with a maximum 
change between two (large) time steps being less than 5*10 ^ . 


CONCLUSIONS 


In this report a computer code is described that can be applied tc the solution 
of an approximate form of the time-dependent Navier-Stokes equations over airfoil 
sections. The governing equations are more general than the conventional boundary 
layer equations, notably in the Inclusion of streamwise and spanwise diffusion terms, 
although the pressure is still imposed by the external flow, as in conventional 
boundary layer theory. The computer code which treats the governing equations incor- 
porates the split LB I scheme in conjunction with QR operator scheme that permits a 
variety of spatial difference schemes, Including standard second order finite differences, 
exponential type methods and fourth order OCI techniques. In the split LBI scheme, an 
implicit sweep is performed in each spatial coordinate direction. A careful ordering 
of these sweeps permits an uncoupling of the continuity equation from the system in 
the first two implicit sweeps. Thus on the first two sweeps the (tridiagonal) system 
block size is reduced from 3 x 3 to 2 x 2 with a resulting cost savings. On the last 
sweep of each time step all the equations in the system are linearly coupled and 3x3 
blocks must be eliminated. Results of several Cartesian problems employing second order 
finite differences indicated that the proposed method is viable and has application to 
more complex flows. Future efforts will aim at exercising the fourth order OCI schemes 
option, treating actual airfoil sections and considering turbulent flow. Comparisons 
will also be made among the various types of spatial schemes to assess the overall 
efficiency gains made by employing higher order methods. The treatment of inflow and 
outflow boundary conditions will also be of primary interest. 


APPENDIX A 


Linearization Technique 

A number of techniques have been used for implicit solution of the 
following first-order nonlinear scalar equation in one dependent variable 

dcfi/ d\ -F[<p) (Al) 

Special cases of Eq. (Al) include the conservation form if F(<}0 = 1, and 
quasi-linear flow if G(<j>) = 4 ,. Previous implicit methods for Eq . (Al) 
which employ nonlinear difference equations and also methods based on two- 
step predictor-corrector schemes are discussed by Ames (Ref. 41, p. 82) and 
von Rosenburg (Ref. 42), p. 56). One such method is to difference nonlinear 
terms directly at the implicit time level to obtain nonlinear implicit 
difference equations; these are then solved iteratively by a procedure such 
as Newton's method. Although otherwise attractive, there may be difficulty 
with convergence in the iterative solution of the nonlinear difference 
equations, and some efficiency is sacrificed by the need for iteration. An 
Implicit predictor-corrector technique has been devised by Douglas and Jones 
(Ref. 43) which is applicable to the quasilinear case (G = <j>) of Eq . (Al) . 

The first step of their procedure is to linearize the equation by evaluating 
the nonlinear coefficient as F(<J> n ) and to predict values of 4> n+1//2 using either 
the backward difference or the Crank-Nicolson scheme. Values for <j> are 
then computed in a similar manner using F(4 n+ ^^) and the Crank-Nicolson scheme. 
Gourlay and Morris (Ref. 44) have also proposed implicit predictor-corrector 
techniques which can be applied to Eq. (Al) . In the conservative case (F = 1), 
their technique is to define G(4>) by the relation G(4) = q>G ( <4>) when such a 
definition exists, and to evaluate G( 4 > n+ ‘*') using values for computed by 

an explicit predictor scheme. With G thereby known at the implicit time level, 
the equation can be treated as linear and corrected values of 4 > n ** are computed 
by the Crank-Nicolson scheme. 

A technique is described here for deriving linear implicit difference 
approximations for nonlinear differential equations. The technique is based 
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on an expansion of nonlinear implicit terms about the solution at the known 
time level, t n , and leads to a one-step, two-level scheme which, being linear 
in unknown (implicit) quantities, can be solved efficiently without iteration. 
This idea was applied by Richtmyer and Morton (Ref. 29 , P» 203) to a scalar 
nonlinear diffusion equation. Here, the technique is developed for problems 
governed by £ nonlinear equations in £ dependent variables which are functions 
of time and space coordinates. The technique will be described for the three- 
dimensional, unsteady equations. 

The solution domain is discretized by grid points having equal spacings 

12 3 12 3 

in the computational coordinates, Ay , Ay and Ay in the y , y*" and y 

directions, respectively, and an arbitrary time step, At. The subscripts i, j, 

12 3 

k and superscript n are grid point Indices associated with y , y , y and t, 

respectively, and thus <!>”., denotes <KyJ, y?, y?, t n ) . It is assumed that 

i> J > K i J R 

the solution is known at the n level, t n , and is desired at the (n+1) level, 
n+1 

t . At the risk of an occasional ambiguity, one or more of the subscripts 

is frequently omitted, so that 4> n is equivalent to . 

1 >3 » R 

The numerical method employed is quite general and is formally derived for 
systems of governing equations which have the following form: 

<)H(<£)/3f = 2>(<£)+S(0) (A2) 

where <j> is a column vector containing £ dependent variables, H and S are 
column vector functions of <f>, and 2 is a column vector whose elements are 
spatial differential operators which may be multidimensional. The generality 
of Eq. (A2) allows the method to be developed concisely and permits various 
extensions and modifications (e.g., noncartesian coordinate systems, turbulence 
models) to be made more or less routinely. It should be emphasized, however, 
that the Jacobian 3H/d<f> must usually be nonsingular if the ADI techniques as 
applied to Eq. (A2) are to be valid. A necessary condition is that each 
dependent variable appear in one or more of the governing equations as a time 
derivative. An exception would occur if for instance, a variable having no 
time derivative also appeared in only one equation, so that this equation could 
be decoupled from the remaining equations and solved a_ posteriori by an alter- 
nate method. 
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The linearized difference approximation is derived from the following 
implicit time-difference replacement of Eq. (A2) : 

(H n+, ~H n )/At =/?[2> (tf)' 1 *') + s n + , ]+(l -/3)[2! {^ n ) + s n ] < A3 > 

where, for example, H n+1 s H(<J> n+1 ). The form of 5) and the spatial differ- 
encing are as yet unspecified. A parameter 0(0 5 0 £ 1) has been introduced 
so as to permit a variable centering of the scheme in time. Equation (A3) 
produces a backward difference formulation for 0=1 and a Crank-Nicolson 
formulation for 3 = 1/2. 

The linearization is performed by a two-step process of expansion about 
the known time level t n and subsequent approximation of the quantity 
(3(j>/3t) n At, which arises from chain rule differentiation, by (<J> n ^ - 4> n ) . The 
result is 

H n+, = H n +(aH/^) n (^ n+ '-^") + 0(AI) 2 <A4a) 

s n+ ' = s" + tas/^) n W n+l -^")+o(AI) 2 (A4b> 

3(r +l ) = ^.(^. n )+(d + (A4c) 


The matrices 3H/3<f> and 3S/3<J> are standard Jacobians whose elements are defined, 
for example, by (3H/3<J>)^ r B 3H^/ 94> r - The operator elements of the matrix 
3^/3<f> are similarly ordered, i.e., (3.2> /3ij>)^ r = ; however, the 

intended meaning of the operator elements requires some clarification. For 
the q*"* 1 row, the operation (3^>^/3^) n (4i n+ ^ - <p n ) is understood to mean that 
{3/3t^ [«ji(x,y,z,t) 3 } n At is computed and that all occurrences of (3<j> r /3t) n 
arising from chain rule differentiation are replaced by (<j>" +1 - <p^)/At. 

After linearization as in Eqs. (A4) , Eq. (A3) becomes the following linear 
implicit time-differenced scheme: 


(dH n /<3<W n + , ~tf> n ) /At =2>(<£ n ) +S n + ft {d!2> /d<p + dS n /d<t>)(4> n +l-<p n ) (A5) 


n+1 

Although H is linearized to second order in Eq. (A4) , the division by At 
in Eq. (A3) introduces an error term of order At. A technique for maintaining 
formal second-order accuracy in the presence of nonlinear time derivatives is 
discussed by McDonald and Briley (Ref. 8), however, a three-level scheme 
results. Second-order temporal accuracy can also be obtained (for 3 * 1/2) by 

A 

a change in dependent variable to <j> = H(<}>), provided this is convenient, since 
the nonlinear time derivative is then eliminated. The temporal accuracy 
is independent of the spatial accuracy. 

On examination, it can be seen that Eq. (A5) is linear in the quantity 
(ij> - <j> ) and that all other quantities are either known or evaluated at 
the n level. Computationally, it is convenient to solve Eq. (A5) for 
(4> n+ * - 4> n ) rather than This both simplifies Eq. (A5) and reduces 

roundoff errors, since it is presumably better to compute a small 0(At) change 
in an 0(1) quantity than the quantity itself. To simplify the notation, a 
new dependent variable ifi defined by 

s (A6) 

n*f*l xvHL n n 

is introduced, and thus \p =4 - , and =0. It is also convenient 

to rewrite Eq. (A5) in the following simplified form: 

(A+ At i 7 )'/'"*' = A1 (£ n )+S n ] (A7a) 

where the following symbols have been Introduced to simplify the notation: 

-ftAl(dS n /d4>) (A7b) 

Jts ~ft(d 2>/d<f>) (A7c) 

It is noted that J?( i|/) is a linear transformation and thus J!(0) = 0. Further- 
more if .A<j>) is linear, then J(ty) = -&2>(ip). 
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Spatial differencing of Eq. (A7n) is accomplished simply by replacing 

12 11 

derivative operators such as 3/3y , 3 /3y 3y by corresponding finite 
difference operators, D^, D^, Henceforth, it is assumed that 2b and 2 have 
been discretized in this manner, unless otherwise noted. 

Before proceeding, some general observations seem appropriate. The 
foregoing linearization technique assumes only Taylor expandability, an assump- 
tion already implicit in the use of a finite difference method. The governing 
equations and boundary conditions are addressed directly as a system of coupled 
nonlinear equations which collectively determine the solution. The approach 
thus seems more natural than that of making ad hoc linearization and decoupling 
approximations, as is often done in applying implicit schemes to coupled 
and/or nonlinear partial differential equations. With the present approach, 
it is not necessary to associate each governing equation and boundary condition 
with a particular dependent variable and then to identify various "nonlinear 
coefficients" and "coupling terms" which must then be treated by lagging, 
predictor-corrector techniques, or iteration. The Taylor expansion procedure 
is analogous to that used in the generalized Newton-Raphson or quasi- 
linearization methods for iterative solution of nonlinear systems by expansion 
about a known current guess at the solution (e.g., Bellman & Kalaba, Ref. 45). 
However, the concept of expanding about the previous time level apparently 
had not been employed to produce a noniterative implicit time-dependent scheme 
for coupled equations, wherein nonlinear terms are approximated to a level of 
accuracy commensurate with that of the time differencing. The linearization 
technique also permits the implicit treatment of coupled nonlinear boundary 
conditions, such as stagnation pressure and enthalpy at subsonic inlet 
boundaries, and in practice, this latter feature was found to be crucial to 
the stability of the overall method (Ref. 30). 
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APPENDIX B 


Tensor Notation and General Coordinates 

n i 2 3 n 

In tho space R with coordinates (x , x , x , ,x ) let the position 

vector be denoted by r. A basis can be formed by constructing the covariant 

vectors e, 

1 

-e = JIL 

1 <*x' (Bl) 


The inner (dot) product of these vectors forms the function 
refer to as the metric tensor 


hy 


which we will 


«U = 



(B2) 


The reason for this terminology will be described subsequently. In three- 

dimension g. . has the form 
ij 



' g .l 


g .r 

S5 

g 

g 

g , 


y 2l 

22 

23 


A. 

g 32 

g 33_ 


The metric tensor is symmetric in the indices i and j so that there are 6 
independent components. For an orthogonal coordinate system the off diagonal terms 
g^i lj*j are zero. 

The differential arc length ds is defined as 


( ds ) 2 - dx ' dx~ 



- d*} o r x J ej ej 
ds 2 = g.. dy} dv) 


(B3) 
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Expanding we obtain 

ds z = g (jtfx'o'x 1 + Q zz d* z dn z + g 33 £/x 3 c/x 3 ^ 

+ 2 g ^tfx'tfx 2 + zg^tfx't/x 3 + 2g 23 tfx 2 tfx 3 

For orthogonal curvilinear coordinates the arc length reduces to 

ds z = gu^x'^x 1 + g 22 c/x 2 tfx 2 +g 33 tfx 3 tfx 3 


2 2 2 

so that g^, g 22 » 833 can be identified with the scale factors as h^, h 2 and h 3 
respectively. 

The Kronecker delta function is defined as 


"e »ei = 8 


.1 l = i 

0 l ^ J 


(B6) 


Notes on raising and lowering indices 

(0) multiplying by changes the index but keeps it in same position 

rv i 

Uj 8| = u, 

u» 8/ = u j 


ii 

(b) multiplying by g^ and g 


u V U J 

lowers 

index 

ij i 

u,g J = u J 

raises 

index 


(c) multiplying 6^ by the metric tensor g^ 


J J k - 

Sj g - b 

R l J k - n lk 

Sj g J = g 


raises index 
substitutes index 


51 


ORIGINAL PAGE If 
OF POOR QUALITY 


Hence we obtain the following relationship 


a 


Ik 


s 9 


Ik 


Ik 

Note that 6 is not the Kronecker delta function in general. Only in the case of 

lk ik 

Cartesian coordinates where g * 1 and i « k Is <5 the Kronecker delta function. 
In order to obtain the contravarlant basis e we write 


Cj • ej = g,j = g l>( 8j 


Employing equation (B-6) 

■«l •e'j 5 S|k k '®"j ] 


we obtain 

°i = g |k -e k 


(B7) 


Hence the metric tensor 
base vectors. 

We may now define g*^ 


relates the covariant to the contravarlant 
as 


9 


)k 



We wish to relate g*^ to g 


jk‘ 


Dotting equation (B~7) by e J we obtain 


- -J — k 
ej • e = g||< e • e 



<B8) 


(B9) 


Now 


from the definition of g 1 ^ (Equation B-8) we obtain the desired result 


-i -I ij Ik * j _ Jk r J r 

e e = g =g 8 k =ge-e k 
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AK. 

Therefore, g also relates the covariant and contravarlant basis vectors, l,e. 


"e * = g^ej^ (bio) 


ij 

The metric tensors g and g, . can also be nsed to raise and lower indices, 

■LJ Ik «k 1 ii 

For example, consider the product of the tensors Cj k , C J C 4 , and with g J 




or g 




c Jk9 


‘J = -J 


c jk 0|j 


•k 

■* 


c/g'l = c kl 


c k 9ij = o 


Ik 


Derivatives 


The derivative of a basis vector e. is a vector in R so that it can be 
represented as a linear combination of other basis vectors, e. ' s. Let e. be a 

. K i 

vector, then its derivative with respect to x J is 


del _ -k^. 

7^ = D i e l = r jl e k 


^9 j 

Dotting this expression with e we obtain 


D j e i • e " = r ji e k •' = r* s J k = r 'f 


.k— 


where rji is the Christoffel symbol. 

r 

Ji " 


V ^ ~ r\ A 

1 « . = D . e , • e 
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£ 

Note that the Chrlstoffel symbol Tjl Is not a third order teneor. 
Representing Tj! in terras of the g^’s we obtain 

®mj* V T J 

°i = = 0,^.1, + VD,T, 

°i g mj = ^lm g kj + g mk*lj 


Similarly we can show that 


D j 9 m i" 9k | r jm + 9 m k Tj | 


and 



-nk 
4 ' mj 


Employing the symmetry property of the T's in their lower index, i.e. 



and adding the three derivatives we obtain 




,km 


[ D i 9 m i + D, Smi- D m9ii] 


which relates the Chris toff el symbol to the derivatives of the components of 
metric tensor. 
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The Kronecker delta function is defined as 



— k - 
e • ej 


Taking its derivative, we obtain 

D| S| k = 0 = D,"e k • Tj 

o = Dj"e k *Tj +"e k • Dj e j 

Substituting in the expression for the derivative of a basis vector 

■— r-,m — 

Dj ej = ljj e m 


and rearranging, 


we obtain 




• r ij e m 






• e 


j 


which relates the Christoffel symbol to the derivative of a contravariant basis 
vector. 

Other properties of Christoffel symbols are 


(a) 


U 


,k a 


% « ■ r if ■ r ii s ; 


(b) 


r ij \i " ^ij£ 


There is also symmetry of lower indices, i.e. 

r ijk = r jik 
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Consider the derivative of a vector u, expressed in a covariant basis 

"u = U* e’| 

-|j*j = D j u ^ e*j = Dj u'"?! + u f Dje] 

Now the derivative of a basis vector is 

Dje t = Ijj e k 

so that the last term in the expression becomes 

u'djTi = u 1 Ifj k T k 

Since k and i are dummy indices, we may interchange Indices, i.e. 

i — i r k_ k r i - 
u D j e ! = u Ijj e k = u r kj e; 

Combining, we obtain the desired expression 

d u / I k pi \ ♦ 

j- = (Dju + u F k j) ej 


The term in the parentheses is called the covariant derivative and is denoted 


as follows 



In a Cartesian coordinate system the covariant derivative reduces to the partial 
derivative, i.e. 

Consider the derivative of a vector in a contravariant basis which is given by 

D j u.e’ J = D, u.T J + UjD.e' J 
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The last term can be expressed in terms of Christoffel symbols, i.e. 

_J pi— k P k_j 

Uj D| e J = uj Tj, e = u k l^j e 1 

Thus we obtain the desired result 

D,rMi> lV = Cu j(| -u k r u * , )e ) 

Djuje^ = uj|,e j 

Now in Cartesian coordinates the covariant derivative of the metric tensor is zero, 



Since this is a tensor equation it is valid not only in a Cartesian system but in 
every coordinate system, i.e. 


Consider now the vector u^. 


Taking the covariant derivative of u and using the previous relationship, we obtain 

X* 

the desired result 

u lll = <“ k 9/k> |, = uk |i «^k + uk 9j'| 

n b 


9|j|l 


= 0 


u l - u 9 ^k 
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An alternate expression for the derivative of a covariant vector is given by 

D i uj'e ) = Uj ||"e = u 1 1| 

-- ( 1 * 1 , 7 , 


Note that the dummy indices are interchanged. 

Vector Operators 

In this section we consider the vector operators, gradient and divergence 
operating on scalars, vectors and tensors. 

Gradient 

Scalar (<£) 

The gradient of a scalar is given as 

Vu =T k ®D k c£ = g km <T m ® D k <£ 


where ® denotes a tensor product. 
Vector (7 = JT| ) 

The gradient of a vector becomes 


Vu = "e k ® D k u e, = e k ®{D k u i e, + u 1 D k e ( } 


= u 1 j k "e” k ® ej 
= g mk u'| k Tm®ei 
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which is a second order tensor. 

Tensor 

Consider a tensor of the form T» = ® ej . The gradient of this tensor 

can be obtained as follows 

VU = e^® ej ® ej 


="e k ® {D k a ij (T, ®ej ) + a lj, (D k e, )®ej +0'^ ®D k ej} 
="e k ®{D k a^(ej ®e[ ) + ©ej + a^lij ® IJJj e^} 


which yields a third order tensor 

| 

i 

Vff = {D k a | i + 0 ( ir K l i + a l4 r^ e k ®ej®e ' 

t 

= { v" + a * 1 + a 14 ^} g mk e~ ® 7 , ® eT j 

= ( a IJ | k ) g mk e m ®e;®ef | 

H 

j; 

if. 
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Divergence 


The divergence of a vector of the form (u = U* ej ) can be obtained as folic -s 


V-TT » e k • D k U l ej 

="e k » u 1 | k “e - ! a u* | k "e k *e*| 

= u'| k Sj k 
V*TT =u k |. 


From the definition of the covariant derivative we obtain 

uVv k + u 'r t k k 


It can be shown (Ref. 13) that 



Substituting this expression into the definition of the covariant derivative we 
obtain an alternate form of the divergence, involving the Jacobian 


Ik 8 D k u + ^ T D i J 


V 


+ u 


k I 


V 


TV Jukl 


Note that this relationship is in conservation form. 
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Consider now the tensor of the form {U = a lj ?, ® ) 

Again we obtain Its divergence as follows j 


V--U = V-a'i^i ®ej 

s e* k ' ( a ^ e“j ® ) 


= e k [D k a lj e) ®e] + a lj D R e", ©e] + a ij T,®D k ej] 
B D k° IJ e“ k * ej® e] + a iJ ^e k ^ ® e] + a iJ e k • ej® T k j 




+ ° ij n.f 8? T. + a ij r* 8 k 


ki u i e j 


by interchanging indices j and l we obtain 


V • u 



+ a 


lj pi 
l l\ 



which reduces to 
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or 

V’ u = { j) + a U r,i]e] 


We will consider now the tensor of the following form ( u - a| e j ®1T^) 
The divergence of u becomes 


V- u 



® e4 aj D k ej ®e ^ + Qj e*! ®D k e^] 


= e 



® el+a ! r kNi® T) - a j e i® r u' i '] 


Rearranging and interchanging dummy indices we obtain 

V. u = [ D h o' J + aj r k ' 4 - a} r' ] (e k • e, ) ® e 1 

= M ■*-afr l l ,-aUf]e ) 


Now the contravariant basis vector is defined as 


e 


i 



Q 


jk 


Substituting into the expression for V*u and with some manipulation we obtain 
the divergence of the second order tensor as 


V • = [ Dj a| - oi ] g \ 

v -“- 
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In order to obtain V*u in the coordinate system we employ the property of the 
metric tensor to raise and lower Indices 

v '“ “ [ D k°| + “f r H" °!t ^hj ®«i 




and obtain the desired result, 


v,u = {[ D k°i 
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TABLE I. - OPERATOR COEFFICIENTS FOR STANDARD 
OPERATOR COMPACT IMPLICIT SCHEME 



where 


P j " hb j 
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TABLE II. - OPERATOR COEFFICIENTS FOR GENERALIZED 
OPERATOR COMPACT IMPLICIT SCHEME 

qf ■ 6 + [p,-3]RCj + [ P 2 ]rcj* 
q,‘ - 60 + [ I0 P |] Rcj + [ P 3 ] nc i + [ T J*|P.] Re * 
q/ ■ G +[p,+3]rcj + [p, + pJ R 9j* + [pJ R t) S 

vhere 

p, • 3 , p B - O , p s ■ mox[ir |t Tr ? ] 


it, • (Tj.j + T^Jpj + r J+|Pl + tt, ir t . l5-2p I +(o- z -l)p ( -3(T,, | + <r !! ) + ir l 


7 r. 


0 

o 

A3 

b~ 

7T_ ■ 

0 

2p,“or 2 >0 

-|" <T| Z (10 - T J+| - T j_, ) 

o-,<°] 

Z 

(2p r o 2 ) z /B 

2p,-cr l <0 


2 


, (T H~V 

,0 "‘ r |+i'' T j-i 

2cr z - 3 tj 4| - Tj_j + 10 + 2hrj_, 

l + T ).i] 

7T, ■ Pj-^ + W, +2v J .,(2 + h- 


' b j-. ' 


fH. 

b J-. 


with b sufficiently small so that 


I0bj-bj_ ( -bj 4| > 0 ond 2 + hCj 4i /bj 4| > 0 for J«2,---,Jand Cj^O 

where 

T J-, “ b J-l /b J » T J*‘ ' b J+l /b j 0nd Rc J * hb J 

r j’ r j’ T J given * T1 TABLE * 
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TABLE III.- OPERATOR COEFFICIENTS FOR 

ALLEN SOUTHWELL EXPONENTIAL SCHEME 


r i" 

= Rc. e Rci / ( 1 - e Rc ^) 

'l 4 

* Rcj /(l-e" Rc i) 

f l c 

= — Rcj + cj 

q f 

= 0 

q i c 

= i 

< 

- 0 


where Rc. = hbj 






TABLE V. - OPERATOR COEFFICIENTS FOR SECOND ORDER 

FINITE DIFFERENCES WITH ARTIFICIAL VISCOSITY 



where S = I for RCj > 2 
S = - I for Rcj < -2 
t - -2- 

|RCj| 

and RCj = hbj 

for |Re.| <2, r', r=, r+ qj. qj , Qj 
reduce to standard finite differences. 
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m 'fv w zt. 





w 
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for Oscillating Flat Plate. 
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